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ABSTRACT 



High-resolution N-body simulations of four popular Cold Dark Matter cosmologies (LCDM, 
OCDM, QCDM, and tilted SCDM), each containing 10^ clusters of galaxies of mass Ml. 5 > 
5 X lQ^^h~^MQ in a Gpc^ volume, are used to determine the evolution of the cluster mass function 
from z = 3 to z = 0. The large volume and high resolution of these simulations allow an accurate 
measure of the evolution of cosmologically important (but rare) massive clusters at high redshift. 
The simulated mass function is presented for cluster masses within several radii typically used 
observationally (i? — 0.5, 1.0, and 1.5 /i~lMpc, both comoving and physical) in order to enable 
direct comparison with current and future observations. The simulated evolution is compared 
with current observations of massive clusters at redshifts 0.3 < z < 0.8. The rim=l tilted SCDM 
model, which exhibits very rapid evolution of the cluster abundance, produces too few clusters 
at z > 0.3 and no massive clusters at z > 0.5, in stark contradiction with observations. The 
r2m=0.3 models — LCDM, OCDM, and QCDM — all exhibit considerably weaker evolution and 
are consistent with current data. Among these low density models, OCDM evolves the least. 
These trends are enhanced at high redshift and can be used to discriminate between flat and 
open low density models. 

The simulated mass functions are compared with the Press-Schechter approximation. Stan- 
dard Press-Schechter predicts too many low mass clusters at z = 0, and too few clusters at 
higher redshift. We modify the approximation by a simple parameterization of the density con- 
trast threshold for collapse, which has a redshift dependence. This modified Press-Schechter 
approximation provides a good fit to the simulated mass functions. 

Subject headings: galaxies:clusters:general - cosmology:theory - dark matter - large-scale structure of 
universe 



The local abundance of clusters of galaxies 
places a powerful constraint on cosmological pa- 
rameters: as^m^ — 0.5, where cts is the rms mass 
fluctuation on an 8ft.~^Mpc scale and fim is the 
present cosmological density parameter (Henry & 
Arnaud (1991), BahcaU & Cen (1992), White, Ef- 
stathiou, & Frenk (1993), Eke, Cole, & Frenk 

(1996) , Viana & Liddle (1996), Kitayama & Suto 

(1997) , Pen (1998)). This constraint, while pow- 
erful, is degenerate in and erg. The evo- 
lution of cluster abundance with redshift, espe- 
cially for massive clusters, breaks this degeneracy 
(e.g., Peebles, Daly, & Juszkiewicz (1989), Ouk- 
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Introduction 



bir & Blanchard (1992, 1997), Eke, Cole, & Frenk 
(1996), Viana & Liddle (1996), Bahcall, Fan, & 
Cen (1997), Fan, BahcaU, & Cen (1997), Carl- 
berg, Yee, & EUingson (1997), Henry (1997, 2000), 
Bahcall & Fan (1998), Eke, et al. (1998), Donahue 
& Voit (1999)). The evolution of high-mass clus- 
ters is strong in = 1, low-cg Gaussian models, 
but is much weaker in low-density, ~1 mod- 
els. Therefore, the evolution of cluster abundance 
provides one of the most powerful methods to de- 
termine both and cg- Various analyses us- 
ing this method have shown that the relatively 
high abundance of observed clusters to z < 0.8 
indicates that the mass-density of the universe is 
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low: ~ 0.2 - 0.5, and as ^ I (sec references 
above). However, the total number of clusters 
currently studied at high redshifts is still small, 
and the model comparisons arc generally based on 
the Press-Schechter (1974) approximation which, 
while surprisingly accurate, is known to have bi- 
ases. Direct comparisons with simulations have 
also been used, but these are generally too small to 
find the rare high-mass clusters, especially at large 
redshifts. New observations are currently under- 
way to determine more precisely the evolution of 
cluster abundance using optical. X-ray, Sunyaev- 
Zel'dovich, and gravitational Icnsing surveys. At 
the same time, larger and more accurate cosmolog- 
ical simulations are needed in order to accurately 
determine the expected mass function of clusters 
of galaxies and its evolution with time, in order to 
allow for proper comparison with the upcoming 
observations. 

In this paper we present gigaparsec-scale sim- 
ulations of four popular CDM models which al- 
low an accurate determination of the evolution of 
the cluster mass function from redshift z = 3 to 
z = 0, including the rare and cosmologically pow- 
erful massive clusters. The rapid increase of avail- 
able computing power and new algorithms which 
run efficiently on parallel machines make possi- 
ble a combination of large simulation volume and 
high resolution. For comparison, the models pre- 
sented here are of a larger volume than the Gover- 
nato et al. (1999) simulations, and attain consid- 
erably higher resolution than the Virgo "Hubble 
Volume" simulations (Evrard 1998; Jenkins et al. 
2000). Furthermore, we are able to cover several 
different models, including the first large volume, 
high resolution simulation of tilted Clm = 1 as well 
as a quintessence model. 

The aim of this paper is to provide the foun- 
dation for comparisons with the upcoming obser- 
vations. Thus the evolution of the mass func- 
tion is presented in a manner that can be directly 
compared with observations, using cluster masses 
within a variety of specific radii typically used in 
observations, rather than using the not so easily 
observable virial mass. We discuss the simulations 
and the method of cluster selection in §2. The evo- 
lution of the resulting cluster mass function is pre- 
sented in §3. Comparison with the Press-Schechter 
approximation is made in §4, and effects of resolu- 
tion are discussed in §5. Comparisons with current 



observations are shown in §6, and conclusions are 
summarized in §7. 

2. Creating the Simulated Mass Function 
2.1. The N-body Simulations 

Four currently favored variants of the CDM 
model were chosen, with parameters consistent 
with numerous observational constraints (Bahcall 
ct al. 1999). These observations tend to favor a 
low density universe with fim — 0.3: in addition to 
spatially fiat (LCDM) and open (OCDM) models 
wc include a quintessence model (QCDM). Like 
LCDM, the QCDM is made spatially flat by in- 
cluding a component with negative pressure, but 
the Q component is dynamically evolving and spa- 
tially inhomogeneous (Caldwell, Dave, & Stein- 
hardt 1998). The equation of state of the Q com- 
ponent, w = —2/3, was chosen to be in the range 
favored by observations (Wang et al. 2000). A 
standard Qm = 1 CDM model was also run, with a 
strongly tilted power spectrum so as to avoid over- 
production of present-day clusters. All models are 
normalized (by as) to match the observed present- 
day cluster abundance and be consistent with the 
COBE microwave background normalization. The 
cosmological parameters for the different runs are 
listed in Table 1. 

The transfer function for a given model was 
computed using the CMBFAST code (Seljak & 
Zaldarriaga 1996), modified to handle a dynamical 
energy component in the QCDM case (Caldwell, 
Dave, & Steinhardt 1998). The resulting power 
spectrum was used to generate initial conditions 
by perturbing particles on a rectilinear grid; this 
step was performed with the COSMICS software 
package (Ma & Bertschinger 1995; Bertschinger 
1995). Each simulation was begun at a redshift 
when the RMS density fluctuation was 15%. For 
QCDM, the fluctuations in the Q component are 
negligible by the time the simulation begins (z « 
30). Thus, once the matter power spectrum has 
been computed, the Q component has an effect 
only through the overall expansion factor a{t). 

All four runs contained 512^ = 134 million par- 
ticles in a periodic cube. For the Om = 0.3 runs 
the box length is 1000/i"^Mpc; in the SCDM run 
it is smaller, such that the volume is reduced by 
a factor of 0.3 relative to the other models (see 
Table 1). With this choice of parameters all mod- 
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els, including SCDM, have an individual particle 
mass of 6.2 x 10^^h~^MQ. Also, the spline ker- 
nel softening length in all cases is e = 27/z~^kpc; 
thus there is the same mass and force resolu- 
tion in all the models, simplifying comparisons. 
The spatial resolution is comfortably smaller than 
the ^ 100ft-~^kpc characteristic core size of clus- 
ters, and the high mass resolution assures that 
two body relaxation is unimportant in the clus- 
ter cores. 

The N-body evolution was carried out with the 
TPM code (Bode, Ostriker, k Xu 2000). The sim- 
ulations ran on a 128 processor SGI Origin 2000 
at NCSA (as part of the Grand Challenge Com- 
putational Cosmology Partnership); a single run 
took 11 to 14 hours to complete. A mesh with 
512 cells on a side was used for long-range inter- 
actions, while interactions internal to dense peaks 
were computed with a tree code. The time step 
for the Particle-Mesh portion of the code is set 
such that in a given step 1) the cosmological ex- 
pansion factor a changes by less than 1%, and 2) 
no particle being evolved solely by the PM code 
moves more than 1/10 the cell size. Each tree 
has its own, possibly smaller, time step. This 
was set such that, for at least 97.5% of the parti- 
cles in the tree, the time step was shorter than 
0.05 * MAX(pO-33,e)/MAX(w,w™,), where pc is 
the density of the cell containing the particle and 
Vrms is the rms velocity of the entire simulation. 
The LCDM run took 469 PM steps from .z = 24 to 
2: = 0, and the most massive tree (which contains 
one of the denser peaks) took 1887 steps. 

In the TPM code, only cells above a given den- 
sity threshold pthr are treated at full resolution, 
so care must be taken to deal only with objects 
which lie above this limit. The density threshold 
parameter was set to pthr = 0.9/j + 4.0(t, where a 
is the dispersion of the cell densities; since cr rises 
over time, pthr is initially only slightly larger than 
p and becomes larger as structure forms (see Bode, 
Ostriker, & Xu (2000)). The evolution of pthr is 
similar in the three riin=0.3 models, rising from 
3p at z=10 to lOp at z—1; at z—0, pthr = 17p for 
LCDM and slightly less for the other two mod- 
els. Since structure formation happens later in 
the SCDM model, pthr is much lower, rising from 
4p at z=l to lip at z=0. In this paper we con- 
sider collapsed objects above 5 x W^-^hr^MQ, or 
80 particles within a radius smaller than the code 



cell size, which implies a density well above pthr at 
all times. In the LCDM run, the smallest object 
picked out for full resolution treatment contained 
28 particles at z=0.5 and 42 particles at z = 0; 
the other runs located even smaller objects. Thus 
we are confident that resolution is not a problem 
for the objects we are considering. This will be 
discussed further in §5. 

The present models are the largest high- 
resolution N-body simulations of OCDM, QCDM, 
and tilted SCDM models currently available. 
Compared to the Virgo "Hubble Volume" simula- 
tions of LCDM and rCDM (Evrard 1998; Jenkins 
et al. 2000), the simulated box size Lbox of these 
models is a third smaller, but mass and force res- 
olution is superior by a factor of three. Governato 
et al. (1999) have carried out an OCDM model 
with cosmological parameters similar to the run 
presented here and mass resolution a factor of 
three higher, but with Lbox = 500/i~^Mpc. In §5 
we will discuss one additional higher resolution 
run of the LCDM model, which uses more than 
10^ particles. 

2.2. Cluster Selection 

Once a simulation is complete, it is necessary 
to identify collapsed objects; this section describes 
the algorithm used to identify clusters. For a given 
cluster we arc interested in locating the total mass 
within a set radius Rc ~ l/i~^Mpc, either comov- 
ing or physical. 

We begin by creating a list of possible clus- 
ter centers, using the HOP package developed by 
Eisenstein & Hut (1998). This code calculates a 
density for each particle, and associates particles 
with nearby neighbors possessing a higher density; 
thus a particle at a density maximum becomes the 
nucleus of a group. HOP tends to break up larger 
clusters into a number of smaller groups; we do 
not attempt to reunite these with the REGROUP 
portion of the HOP package but instead all group 
positions are treated as potential cluster centers. 
However, groups with less than 25 particles or with 
central density pc < 180 are eliminated, because 
including them in the analysis only results in more 
poor clusters with mass less than 2 x 10^^h~^MQ, 
well below our mass threshold. 

Nearby pairs of clusters are merged together, 
working out to a separation of lh~^ Mpc; the posi- 
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tion of a merged cluster is set to that of the denser 
of the two progenitors. After this merging is com- 
plete, the position of each cluster is revised to the 
center of mass of all particles located within lh~^ 
Mpc of the original position. Using this final set 
of cluster positions, all particles within a specified 
radius Rc of a cluster center are assigned to that 
cluster. If a particle is within Rc of multiple clus- 
ters, it is assigned to the group to which it has a 
greater binding energy (using the previously com- 
puted mass and central position of each group). 

Note that these cluster masses, defined specifi- 
cally within a fixed radius Rc (for example R,. = 
1.5/i~^Mpc) are not the same as virial cluster 
masses. For a rich cluster with Ryir > 1.5/i~^Mpc 
(where Ryir is the radius containing a mean den- 
sity of 178 times the critical value). Mi. 5 will be 
smaller than the virial mass and for poor sys- 
tems A'/i.5 will be larger than the virial mass when 
Rvir < 1.5/i~^Mpc. Using a fixed Rc facilitates 
a proper comparison with observations, since the 
virial radius cannot easily be determined observa- 
tionally. 

3. Evolution of the Cluster Mass Function 

Using the selection criteria of §2.2, the cluster 
mass function (MP), which represents the num- 
ber density of clusters above a given mass thresh- 
old, is determined from the simulations as a func- 
tion of redshift for z = 0, 0.17, 0.5, 1, 2, and 
3. At z = 0, approximately 10^ clusters with 
M1.5 > 5 X lO^^/i~^M0 are located within the sim- 
ulation volume. The MP is determined for cluster 
masses within different radii — Rc = 0.5, 1, and 
1.5/i~^Mpc, both comoving and physical — in or- 
der to explore the dependence of the MP on this 
scale, as well as provide for proper comparisons 
with different observations. 

The cluster MP is presented in Figure 1 for the 
four models studied. The cluster mass used here 
is the mass within the 'standard' 1.5/i~^Mpc co- 
moving radius, frequently used in observations. At 
z = 0, all models yield the same MP, designed to 
be consistent with observations; this consistency is 
essentially forced by the selected erg normalization 
of each model (Bahcall and Cen 1992, White et al. 
1993, Eke et al. 1996, Pen 1998). The models are 
thus degenerate (by construction) at z r^O, at least 
for the richer clusters (SCDM shows many more 



poor clusters). At higher rcdshifts, a negative evo- 
lution of the cluster MP is seen in all models, as 
expected, reflecting the growth of clusters with 
time from small density fluctuations (e.g., Press 
& Schechter 1974, Peebles 1993, Eke et al. 1996, 
Fan et al. 1997). 

The rate of the evolution, however, is strongly 
model dependent; it breaks the degeneracy seen at 
z ~ 0. The ^^=1 tilted SCDM model, with much 
stronger evolution, predicts orders-of-magnitude 
fewer clusters at high redshift than do the low- 
density models. For example, for a given mass 
the abundance of clusters in SCDM at z = 0.5 is 
comparable to the abundance at z = 1 in the low- 
density models. Similarly, the number of z = 1 
clusters in SCDM is comparable to the number 
of z = 2 clusters in the low-density models. At 
z = 0.5, the abundance of Coma-like clusters (with 
Mi.scom =i 6 X lO^/i-^Mg) is 100 times larger 
in the low-density models than in SCDM. This 
large effect, occurring at relatively low redshifts of 
z ~0.5 1, provides a powerful method for deter- 
mining and as from the observed abundance of 
high-redshift massive clusters. The dependence of 
the cluster abundance on A is much weaker; it be- 
comes more significant at higher redshifts (z >1) 
(e.g. Eke et al. 1996, 1999, Bahcall et al. 1997, 
Fan et al. 1997). 

The LCDM and QCDM mass functions evolve 
at roughly the same rate, consistent with the pre- 
diction of Wang & Steinhardt (1998). QCDM 
evolves somewhat slower; the differences become 
more pronounced at z > 1. This can be quan- 
tified by considering the ratio of cluster abun- 
dance n{z = 1)1 n{z = 0) for clusters of mass 
> 8 X IQ'^^h^'^MQ (using the best- fit P-S approx- 
imation of §4). QCDM with a ratio of 2.2 x lO"^ 
is roughly double LCDM (1.2 x lO"''') and QCDM 
(1.1 X 10^3) (SCDM is much lower at 6.9 x 10"^; 
sec §6). With improved data sets, this will enable 
a discrimination between flat and QCDM cosmolo- 
gies (assuming is known) . The evolution of the 
cluster MP presented in Figure 1 can be used for 
direct comparisons with observations. 

The analysis presented above is for clus- 
ter masses within the standard Abell radius of 
1.5ft.~'^Mpc comoving. Many observations, but 
not all, determine cluster masses within this radius 
(based on velocity dispersion, gas temperature, or 
gravitational lensing). Some observations however 
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determine cluster masses within other radii (typi- 
cahy 0.5 or l/i~^Mpc). In order to provide proper 
comparisons with observations, as well as under- 
stand the dependence of the MF on the selected 
radius within which the cluster mass is measured, 
we investigate the MF evolution for different radii. 
In Figure 2 we present the evolution of the cluster 
MF for cluster masses measured within radii of 
0.5, 1, and 1.5/i~^Mpc, both comoving and physi- 
cal, for each of the four models. 

The dependence of the MF evolution on ra- 
dius is evident in Figure 2. The progression from 
large to small radii for either the comoving or 
physical scales reflects mostly a progression from 
the large cluster mass within 1.5/i~^Mpc to the 
smaller mass contained within 0.5h~^M.pc for the 
same clusters. Therefore, this progression simply 
shifts the MF to smaller masses as the radius de- 
creases (at all redshifts). If the cluster density pro- 
file is p{R) ~ on these scales, where 8—7 « 1, 
then the mass inside R is M{< R) ^ R^~'^ and a 
similar shift of mass with radius is expected — as 
is indeed seen in the simulations. 

A more interesting trend is seen between co- 
moving and physical radii. The evolution of the 
MF is much weaker for cluster masses determined 
within physical radii than for comoving radii. The 
strong evolution seen for comoving radii is greatly 
reduced for physical radii because the same mass 
is contained within the larger physical radius (at 
high redshift) instead of the smaller comoving ra- 
dius (since Rphy = Rcom(l + z)). Therefore, for a 
fixed cluster mass (as a function of redshift) the 
use of physical radius represents relatively lower 
mass clusters at high redshifts as compared with 
the comoving radius, by a factor of approximately 
(1 -|- 2;)^"'*'. Since lower mass clusters are more 
numerous, there will be more clusters at high red- 
shifts for cluster masses defined within physical 
radii than within comoving radii, thus yielding 
considerably weaker evolution for the former. The 
effect is stronger in low density models, where 
only mild evolution is seen in physical coordinates. 
This is mostly due to the fact that for a low density 
universe evolution is significant from early times 
until z ~ ri~^ ~ 3, when the universe becomes 
"open" for local observers and both accretion and 
merging are reduced. In such models, clusters viri- 
alize at z ~ 2, after which their evolution slows 
down; in fixed coordinates their properties, includ- 



ing the mass function, approach constancy. 

The difference is large; for example, the 
abundance of Coma-like clusters (Mi. 5 ~ 6 x 
lO^'^h^^ Mq) at z = 1 is nearly 100 times larger 
for R = 1.5/i~^Mpc physical radius than for 
1.5h~^Mpc comoving radius. It is therefore es- 
sential that observed and simulated clusters be 
compared using the same radii. The results pre- 
sented in Figure 2 can be used for comparisons 
with observations of cluster masses within these 
different radii. 

4. Comparison with Press- Schechter Ap- 
proximation 

In this section we compare the simulation re- 
sults with the predictions of the Press & Schechter 
(1974) approximation frequently used (e.g. Eke, 
Cole, & Frenk 1996, Fan, Bahcall, & Cen 1997) 
in the estimation of the number density and mass 
distribution of collapsed objects. The P-S approx- 
imation states that, after smoothing on the appro- 
priate length scale, regions with density contrast 
above some threshold Sc will have collapsed and 
virialized. The spherical collapse model predicts 
Sc = 1.686, with a weak dependence on f2„ and 
A (Eke, Cole, & Frenk 1996; Wang & Steinhardt 
1998). Because P-S theory predicts the mass in- 
side a virial radius, and this radius is often used 
in numerical work, we first compare the standard 
P-S predictions with the simulated mass function 
found using HOP. Employing the HOP regrouping 
algorithm and tracing clusters to an outer over- 
density threshold of 160 gives results similar to 
the Friends-of-Friends (FOF) algorithm with link- 
ing parameter 0.2, which is frequently used to 
analyze numerical simulations (Eisenstein & Hut 
1998; Governato et al. 1999; Jenkins et al. 2000). 
The resulting HOP mass function is shown in Fig- 
ure 3, along with the P-S prediction. The shape 
of the standard P-S mass function does not fit 
the numerical results; it predicts too many low 
mass clusters and too few higher mass clusters. 
This is now a well established drawback of the P- 
S approximation (Gross et al. 1998; Governato et 
al. 1999; Jenkins et al. 2000). Furthermore, at 
higher redshifts there are more collapsed objects 
than predicted by standard P-S. 

To explore the comparison in more detail, and 
to allow for proper comparison with observations. 
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wc will use the mass function defined in terms of 
the observable quantity of mass within a specific 
radius. Because P-S theory predicts the mass in- 
side a virial radius, an adjustment needs to be 
made to predict the mass within 1.5/i~^Mpc. We 
follow the prescription of Carlberg et al. (1997) re- 
lating the top-hat smoothing length to A^i.s rather 
than the virial mass. A slope of 3 — 7 = 0.6 is 
assumed in the relation M(< R) oc R^~^ near 
the virial radius. This slope is consistent with the 
mass profile of observed rich clusters (Carlberg, 
Yee, & Ellingson 1997; Rines et al. 2000). Also, 
since N-body simulations give results in variance 
with P-S, we treat 6c as a free parameter. For 
a given comoving Mi. 5, 5c was adjusted so that 
the P-S prediction for n(> Mi. 5) matched the N- 
body result at that mass. This was done for all 
models and redshifts, beginning at the mass of the 
tenth most massive cluster and working down to 
6 X lOi3/i-iM0. 

The resulting 6c as a function of Mi.scom is 
shown in Figure 4. It is apparent that the stan- 
dard P-S approximation (using 6c from spheri- 
cal overdensity in linear theory) does not fit the 
simulations. The cluster abundance depends on 
5c roughly as exp(— J^/cr^), so a higher 6c means 
fewer clusters have formed. The 6c required by the 
simulation MF varies with mass. In other words, 
the shape of the standard P-S mass function is 
incorrect; for the low flm models it predicts too 
many low mass clusters (since the standard 6c is 
not large enough), even if 6c is fixed to match the 
high mass clusters (> 2 x IO^/i-^Mq). The tilted 
SCDM model shows the opposite trend — we are 
finding more low mass clusters in the simulations 
than is predicted. This is due to the assump- 
tion that Af(< i?) (X i?"-^ in adjusting the virial 
mass to Ml. 5. For = 1 the extrapolation from 
1.5/i~^Mpc to Ryir reaches to radii < l/i~^Mpc 
for low mass clusters. Allowing a steeper slope 
(M(< R) oc R) for smaller SCDM clusters with 
a virial radius less than l/i~^Mpc yields results 
similar to those seen in the other models. 

A second trend seen in Figure 4 is that 6c is 
lower at higher redshifts, showing that there are 
more collapsed objects at 2 > than predicted by 
standard P-S. The redshift dependence of 5c can 
be parameterized as 6c = 60 + 6i/{l + z), with 5o 
and 61 chosen separately for each model. A lin- 
ear fit was made to the best values at Mi.5com = 



2 X lO^'^h'^ Mq, and then the parameters were 
adjusted slightly to reduce the difference between 
the P-S predictions and simulations to below 10% 
across the entire mass range, if possible. The final 
choices of 60 and 5i are shown in Table 2. It can be 
seen that 5o is close to the canonical value of 1.68 
and that 6^ is < 10% of this value, meaning the 
z dependence is weak. Figure 5 shows both the 
N-body MF and the results of the modified P-S 
formula using the 6c from Table 2, along with the 
fractional difference between the two. The simu- 
lation data are fit well by this modified P-S re- 
lation, to within 10% for z < 2 (or z < 1 and 
M > lO^^h-^MQ for SCDM). 

An alternative to the simple power-law density 
profile used here is the proposed universal profile 
for dark matter halos of Navarro, Frenk, & White 
(1997) (see also Lokas & Mamon (2000)); the "toy 
model" of Bullock, et al. (2000) provides the mass 
and redshift dependence of the concentration. In 
order to test whether including this dependence 
would improve the analytic fit, the above analysis 
was repeated using the NFW profile to adjust the 
virial mass to Mi. 5. No improvement in the fit was 
seen; 6c shows both mass and redshift dependence 
in a manner similar to that observed using the 
power-law profile, and the fit is somewhat worse 
for redshifts z '>2. 

Our results are in agreement with Governato et 
al. (1999), who likewise found the best fit 6c to be 
slightly lower as redshift increases. However, the z 
dependence found here is weaker. Governato et al. 
(1999) simulated an open model with parameters 
close to our OCDM run, and found 6c ~ 1-65 at 
z=l, agreeing well with our value. However, at 
z=0 the Governato et al. (1999) value is well above 
ours; this may be due to the fact that their erg 
value is higher as well. 

5. A Higher Resolution Run 

An additional simulation was carried out us- 
ing the same parameters as the LCDM run, ex- 
cept with 1024^ particles, allowing an examina- 
tion of the effect of increasing the numerical res- 
olution. The increased number of particles re- 
duces the mass of particle by a factor of eight to 
7.75 X 10^°/i~^Mq; the softening length was cho- 
sen to be 13.6 hr^kpc. The run took 800 PM steps 
and up to 9000 tree steps. This simulation was 
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carried out on a 256 processor SGI Origin 2000 at 
NCSA and took 9.5 days; the considerable compu- 
tational expense precludes further runs of this size. 
For comparison, the N = 1024^ Virgo "Hubble 
Volume" LCDM simulation (Jenkins et al. 2000; 
Evrard 1998) contained 27 Gpc^, with particle 
mass 2.25 x lO^^/i^^M© and Plummer softening 
length 100 h'^kpc. 

A comparison of the mass functions found in 
the 512^ and 1024^ runs is shown in Figure 6. It 
can be seen that for ^; < 2 the results are similar, 
the main difference being that the higher resolu- 
tion run yields slightly more lower mass clusters: 
at 5 X lO"/i-iM0, there is a 10% difference. For 
typical rich clusters, with mass > 10^^/i~^A'/q, the 
effect is negligible. The modified Prcss-Schechter 
prediction, using 6c as given in Table 2, is also 
shown. At 2; = 0, the simulation and modified P-S 
mass functions arc within 5% of each other down 
to n(> Ml. 5) = 10~'^/i3Mpc"3 (that is, when there 
are more than 100 clusters in the simulation vol- 
ume). At higher redshifts the modified P-S pre- 
diction yields roughly 10% fewer low mass clusters 
than the 1024^ simulation. Since increasing mass 
resolution by a factor of eight and doubling force 
resolution has only a minor effect, we feel confident 
the results presented here will not be changed sig- 
nificantly by increasing resolution. 

6. Comparison with Observations 

The evolution of cluster abundance is compared 
with observations in Figure 7. Here we plot the 

abundance of massive clusters, above a given mass 
threshold, as a function of redshift. We use the 
most distant observed massive clusters since such 
clusters provide the most powerful test of cosmo- 
logical evolution: these massive clusters cannot 
exist at high redshifts in Qm — 1 Gaussian mod- 
els; they were simply not formed yet (e.g., Pee- 
bles 1993, Fan et al. 1997, Bahcall and Fan 1998, 
Donahue et al. 1999). The clusters we use are 
the same as those used by Bahcall & Fan (1998): 
MS 0016-M6, 0451-03, and 1054-03. Those clus- 
ters are at z —0.55 and z =0.83, with masses 
above a threshold of Mi.scom > 8 x lO^^/i^^Af©. 
They have measured velocity dispersion (o"r ^ 
1200 km/s), gas temperature (T > 8 Kev), strong 
S-Z decrements (Carlstrom 1998), and (in 2 of the 
3 clusters) gravitational lensing mass determina- 



tions. Also shown arc clusters at z = 0.38 from 
the temperature function of Henry (2000). The 
temperature has been converted to Mi.5com us- 
ing the observed M-T relation: M{< l/i-^Mpc) = 
A;T(keV) x 10^^/i~^Mo, with the minor extrapola- 
tion to i?i.5com following the observed cluster mass 
profile (Hjorth, Oukbir, & van Kampen (1998), 
Bahcall & Fan 1998, Bahcall & Sette 2000). The 
scatter in this observed relation is ± 25%. 

The model predictions are presented by the 
curves in Figure 7. The shaded horizontal band 
indicates the range where only a single cluster is 
expected to be found in the simulation box; be- 
low this band, where no clusters can be found, we 
extrapolate the n{z) curve with our modified P-S 
approximation (§4) as shown by the lighter curves. 
The comparison with observations is powerful: no 
massive clusters are found in the = 1 simu- 
lations at any redshifts z > 0.2! A few clusters 
are found in the Qm = 0.3 simulations at 2; ~ 0.6 
and 0.8, consistent with observations. The fact 
that such massive clusters exist in the universe at 
2: ~ 0.8 and 0.6 rules out fim = 1 Gaussian mod- 
els at a very high significance level: if flm = 1, 
10~® clusters would be expected in the observed 
survey volume at 2 ~ 0.8 while 1 cluster is ob- 
served; 10~^ clusters would be expected at 2; ~ 
0.6 while 2 arc observed; and 4 x 10~^ clusters 
expected at 2: ~ 0.4 while 2 are observed. Equiva- 
lently, based on the observed density, the simula- 
tion box should contain 10 such massive clusters 
at 2 ~ 0.8 and 10 clusters at 2 ~ 0.6, while the 
flm = 1 simulation contains no clusters; from the 
modified P-S extrapolation one would expect only 
lO"'^ and 10~^ clusters, respectively. In the Qm = 
0.3 simulations, on the other hand, we find 1 clus- 
ter in the simulation box at z = 0.8 and 3 clusters 
at 2 = 0.6. An even better match with the ob- 
served cluster abundance would be achieved with 
a somewhat lower fl„i value of ^0.2 (e.g., Bahcall 
& Fan 1998). We note that the observed cluster 
abundance at 2=0.5-0.8 could be even higher than 
used above if the observational surveys are incom- 
plete; this would further reduce the best fit value 

of Qrn- 

The models with ~ 0.3 (LCDM, QCDM, 
OCDM) provide a good fit to the data. While a 
precise value of flm needs to await a larger and 
improved sample of massive high redshift clusters, 
the current data — the existence of just a few mas- 
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sivc clusters at high-rcdshift — already shows that 
(for Gaussian models) the mass density of the uni- 
verse is low: flm ~ 0.3. These conclusions from 
direct large-scale simulations reinforce previous es- 
timates based on the standard P-S approximation 
(Bahcall et al. 1997, 1998, Carlberg et al. 1997, 
Henry 1997, 1999, Eke et al. 1998, Donahue et al. 
1999, and references therein). 

7. Conclusions 

Gigaparsec scale high-resolution N-body sim- 
ulations are used to determine the cluster mass 
function from 2; = 3 to 2; = for four CDM mod- 
els: LCDM, OCDM, QCDM, and tilted SCDM. 
These large, high-resolution simulations, with 10^ 
clusters of mass Mi. 5 > 5 x IO^^H'^Mq at the 
present day, allow an accurate determination of 
the MF evolution over this redshift range. 

The evolution of the mass function is presented 
for cluster masses within several radii typically 
used in observations: 0.5, 1, and 1.5h~^Mpc, both 
comoving and physical; this will enable direct com- 
parisons with different observations. The evolu- 
tion of the MF for clusters of a given mass within a 
physical radius is found to be considerably weaker 
than for mass within a comoving radius; this is 
due to the relatively larger physical radii at high 
redshift, which means that more numerous poorer 
clusters are included. After virialization at 2; « 2 
in the flm = 0.3 models, there is relatively little 
evolution in physical coordinates. Because of the 
large differences in the resulting evolution, it is 
essential that the cluster evolution be compared 
using the same simulated cluster mass — within 
the same radius — as observed. 

We compare the simulation results with the 
Press-Schechter approximation; the latter predicts 
too many low mass clusters at 2 = and too few 
clusters at higher redshifts. We provide refined 
parameters for this approximation that best fit 
the simulations as a function of redshift and mass. 
With the chosen fits to the ovcrdcnsity parame- 
ter Sc, the agreement with the simulation MF is 
better than 10% for 2 < 2 in the three flrn = 0.3 
models and z < 1 in SCDM. 

We compare the simulation results with current 
observations of the most massive distant clusters 
observed to z ~0.8. We show that Gaussian 0^=1 
models such as tilted SCDM, which predict ex- 



tremely rapid evolution from z ^1 to z ^0, are 
ruled out at a high significance level; they pro- 
duce no massive clusters at high redshifts (10~^ 
clusters per the observed volume), in contradic- 
tion with observations. Showing much weaker evo- 
lution, the flm = 0.3 models— LCDM, OCDM, 
and QCDM — are all consistent with the data. 
These results provide one of the most powerful 
constraints on the mass density of the universe. 
Upcoming surveys of distant clusters should be 
able to provide further discrimination between fiat 
(LCDM or QCDM) and open (OCDM) low den- 
sity cosmologies, since the latter shows relatively 
weaker evolution. 

This research was supported by NSF Grants AST- 
9318185 and AST-9803137 (under Subgrant 99- 
184), and the NCSA Grand Challenge Computa- 
tional Cosmology Partnership under NSF Coop- 
erative Agreement ACI-9619019, PACI Subaward 
766. 
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Mass Functions at z = 0.0, 0.5, 1.0, 2.0 : 




10" 10'= 
M(<1.5 Mpc„J h-' Ms 



Fig. 1. — The cluster mass functions at redshifts 
of 2; = 0, 0.5, 1.0, and 2.0 for the four different cos- 
mological models are plotted as lines. Note that 
the SCDM model (dotted line) a.t z = 0.5 is com- 
parable to the low density models at z = 1 and the 
SCDM model at z = 1 is comparable to the low 
density models at = 2. The three data points 
with error bars are from observational data dis- 
cussed in §6. 

This 2-column preprint was prepared with the AAS MJjX 
macros v5.0. 
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Fig. 2. — The evolution of the cluster MF for clus- 
ter masses measured within radii of 0.5, 1, and 
1.5h~^A/[pc, both comoving and physical, for each 
of the four models. 
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Fig. 3. — Solid lines: the simulated mass func- 
tion, using masses determined via the HOP al- 
gorithm. Dotted lines: Press-Schechter ap- 
proximation, with 5c from linear theory. The 
curves from top to bottom are for redshifts z = 
0,0.17,0.5,1,2,3. 
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Fig. 4. — The best fit P-S overdensity parameter 
6c as a function of mass. The curves from top 
to bottom are for redshifts z = 0, 0.17, 0.5, 1, 2, 3. 
Virial mass is translated to mass within comov- 
ing 1.5/i~^Mpc assuming the observed mass profile 
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Fig. 5. — The larger panels show the number den- 
sity as a function of mass from the simulations 
(solid line) and from the modified P-S approxima- 
tion (dotted lines) using the overdensity parame- 
ter from Tabic 2. The smaller panel beneath each 
of these plots shows the fractional difference be- 
tween the two, n(PS)/n(sim)-l. The line types 
match those in Figure 4. 
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Fig. 7. — The abundance of massive clusters, 
above a given mass threshold, as a function of 
redshift. The dark lines arc from our numerical 
simulations, while the light lines are the modified 
P-S approximations using the parameters from Ta- 
ble 2. The data points are from observations dis- 
cussed in §6. The dotted region indicates the range 
of densities for which the numerical simulations 
would only have one cluster in the entire volume. 



Fig. 6. — A comparison of the mass function in the 
LCDM model using 512^ (dashed hues) and 1024^ 
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Table 1 
Model Parameters 



Model 










hiQo 


n 


OS 




LCDM 


0.30 


0.70 


0.040 


0.260 


0.67 


1.0 


0.90 


1000 


OCDM 


0.30 


0.00 


0.040 


0.260 


0.67 


1.0 


0.80 


1000 


QCDM 


0.30 


0.70^ 


0.040 


0.260 


0.67 


1.0 


0.84 


1000 


SCDM 


1.00 


0.00 


0.076 


0.924 


0.50 


0.625 


0.50 


669.4 



''box length in h ^Mpc 

''actually a Q component with w=-2/3 



Table 2 
Best-fit Sc = So + Si/{1 + z) 



Model 6o Si 

LCDM 1.60 0.18 

OCDM 1.61 0.07 

QCDM 1.61 0.13 

SCDM 1.53 0.09 
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